##########################################################################
## Replication syntax for Malet & Hegewald (2025, JEPP):
## "The Changing Geography of Support for European Integration in the Shadow of the Ukraine War"
## Results manuscript ex mrp
##########################################################################

# R version: 4.5.1 
# Platform: x86_64-apple-darwin20

# here 1.0.1
# readxl 1.4.5
# tidyverse 2.0.0
# reshape2 1.4.4
# eurostat 4.0.0
# fixest 0.12.1
# broom 1.0.9
# ggpubr 0.6.1

# 1. Load packages #######################################################
library(here)
library(readxl)
library(tidyverse)
library(reshape2)
library(eurostat)
library(fixest)
library(broom)
library(ggpubr)

# 2. Basic setup #########################################################
rm(list=ls())
i_am("01_results_manuscript_exmrp.R")

# 3. Import data ##########################################################
data <- read.csv("survey_data.csv")
elec_long <- read.csv("electricity_data.csv")

# 4. Figure 1 ############################################################
mani <- data %>% select(studyno1, UkraineStatements_securCountry, NUTScode, closeness_to_russia) %>%
  na.omit()

mani$UkraineStatements_securCountry <- ifelse(mani$UkraineStatements_securCountry<5, 
                                              mani$UkraineStatements_securCountry, NA)
table(mani$UkraineStatements_securCountry)

mani <- mani %>%
  mutate(UkraineStatements_securCountry = 5 - UkraineStatements_securCountry)

table(mani$UkraineStatements_securCountry)

mani_av <- mani %>%
  mutate(UkraineStatements_securCountry = as.numeric(as.character(UkraineStatements_securCountry))) %>%
  group_by(NUTScode, closeness_to_russia) %>%
  summarise(avg_UkraineStatements = mean(UkraineStatements_securCountry, na.rm = TRUE))

cor_test <- cor.test(mani_av$closeness_to_russia, mani_av$avg_UkraineStatements)
cor_test

cor_coeff <- round(cor_test$estimate, 3)
p_value <- cor_test$p.value

sig_label <- ifelse(p_value < 0.001, "***",
                    ifelse(p_value < 0.01, "**",
                           ifelse(p_value < 0.05, "*", "ns")))


Fig1 <- ggplot(mani_av, aes(x = closeness_to_russia, y = avg_UkraineStatements)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE, color = "blue") +
  labs(
    x = "Logged closeness to Russia \n(rescaled 0-1)",
    y = "Respondent's agreement with statement that \nRussia's invasion of Ukraine is a threat to their country\n ←Disagreement            Agreement→"
  ) +
  annotate("text", x = min(mani_av$closeness_to_russia), y = max(mani_av$avg_UkraineStatements), 
           label = paste0("r = ", cor_coeff, " ", sig_label),
           hjust = 0, size = 5) +  
  theme_minimal()

png(file=here("Figure1.png"),
    width = 6, height = 6, units = 'in', res = 800)
Fig1 
dev.off()

# 5. Figure 2 ############################################################
elec_long$date[elec_long$semester=="2018-S1"] <- "2018-01-01"
elec_long$date[elec_long$semester=="2018-S2"] <- "2018-07-01"
elec_long$date[elec_long$semester=="2019-S1"] <- "2019-01-01"
elec_long$date[elec_long$semester=="2019-S2"] <- "2019-07-01"
elec_long$date[elec_long$semester=="2020-S1"] <- "2020-01-01"
elec_long$date[elec_long$semester=="2020-S2"] <- "2020-07-01"
elec_long$date[elec_long$semester=="2021-S1"] <- "2021-01-01"
elec_long$date[elec_long$semester=="2021-S2"] <- "2021-07-01"
elec_long$date[elec_long$semester=="2022-S1"] <- "2022-01-01"
elec_long$date[elec_long$semester=="2022-S2"] <- "2022-07-01"
elec_long$date[elec_long$semester=="2023-S1"] <- "2023-01-01"
elec_long$date[elec_long$semester=="2023-S2"] <- "2023-07-01"
elec_long$semester <- semester(elec_long$date, with_year = TRUE)

Fig2 <- elec_long %>% 
  group_by(semester, cntry) %>% 
  ggplot(aes(x = semester, y = elec_pr, group = cntry, color = cntry)) +
  geom_line() +
  ylab("Electricity price \n(Euro per kWh)") + xlab("Semester") +
  facet_wrap(~cntry) + 
  theme_minimal(base_size = 14) +
  theme(legend.position="none", axis.text.x = element_text(angle = 90))

png(file=here("Figure2.png"),
    width = 11, height = 8, units = 'in', res = 800)
Fig2 
dev.off()

# 6. Figure 3 ############################################################
data_map <- data %>% select(NUTScode, exposure)

nuts1 <- get_eurostat_geospatial(nuts_level = 1, year=2021)
nuts2 <- get_eurostat_geospatial(nuts_level = 2, year=2021)
nuts3 <- get_eurostat_geospatial(nuts_level = 3, year=2021)
nuts_shapes <- bind_rows(nuts1, nuts2, nuts3)

nuts_shapes <- nuts_shapes %>% rename(NUTScode = id) %>% 
  select(NUTScode, geometry)

data_map <- data_map %>% left_join(nuts_shapes, by = "NUTScode")
data_map <- unique(data_map)

Fig3 <- ggplot(data_map) +
  geom_sf(aes(geometry = geometry, fill = exposure)) +
  scale_fill_viridis_c(
    option = "viridis", 
    name = "Exposure to \nenergy crisis",
    breaks = c(0, 0.5, 1)  
  ) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "bottom")

png(file=here("Figure3.png"),
    width = 12, height = 8, units = 'in', res = 800)
Fig3 
dev.off()

# 7. Figure 4 ############################################################
data_line <- data %>% select(studyno1, date, Proposals_CommonDefence, Proposals_CommonForeign, Proposals_CommonEnergy, Proposals_EnlargeEU) %>%
  mutate(date = ymd(date, quiet = TRUE))

median_dates <- data_line %>%
  group_by(studyno1) %>%
  summarize(
    median_date = median(date, na.rm = TRUE),
    start_date  = min(date, na.rm = TRUE),
    .groups = "drop"
  )

data_line_fin <- left_join(data_line, median_dates, by = "studyno1")

data_line_fin <- data_line_fin %>%
  filter(!median_date %in% as.Date(c("2020-11-07", "2021-09-28"))) 

data_line_fin <- data_line_fin %>%
  group_by(median_date) %>%
  summarise(
    Defence_mean = mean(Proposals_CommonDefence, na.rm = TRUE) * 100,
    Defence_se   = sd(Proposals_CommonDefence,   na.rm = TRUE) / sqrt(sum(!is.na(Proposals_CommonDefence))) * 100,
    Foreign_mean = mean(Proposals_CommonForeign, na.rm = TRUE) * 100,
    Foreign_se   = sd(Proposals_CommonForeign,   na.rm = TRUE) / sqrt(sum(!is.na(Proposals_CommonForeign))) * 100,
    Energy_mean  = mean(Proposals_CommonEnergy,  na.rm = TRUE) * 100,
    Energy_se    = sd(Proposals_CommonEnergy,    na.rm = TRUE) / sqrt(sum(!is.na(Proposals_CommonEnergy))) * 100,
    Enlargement_mean = mean(Proposals_EnlargeEU, na.rm = TRUE) * 100,
    Enlargement_se   = sd(Proposals_EnlargeEU,   na.rm = TRUE) / sqrt(sum(!is.na(Proposals_EnlargeEU))) * 100,
    .groups = "drop"
  ) %>%
  mutate(
    Defence_lower = Defence_mean - 1.96 * Defence_se,
    Defence_upper = Defence_mean + 1.96 * Defence_se,
    Foreign_lower = Foreign_mean - 1.96 * Foreign_se,
    Foreign_upper = Foreign_mean + 1.96 * Foreign_se,
    Energy_lower  = Energy_mean  - 1.96 * Energy_se,
    Energy_upper  = Energy_mean  + 1.96 * Energy_se,
    Enlargement_lower = Enlargement_mean - 1.96 * Enlargement_se,
    Enlargement_upper = Enlargement_mean + 1.96 * Enlargement_se
  )

data_line_fin_long <- data_line_fin %>%
  select(
    median_date,
    Defence_mean, Defence_lower, Defence_upper,
    Foreign_mean, Foreign_lower, Foreign_upper,
    Energy_mean,  Energy_lower,  Energy_upper,
    Enlargement_mean, Enlargement_lower, Enlargement_upper
  ) %>%
  pivot_longer(
    cols = -median_date,
    names_to = c("policy_area", ".value"),
    names_pattern = "(.*)_(mean|lower|upper)"
  ) %>%
  rename(support = mean)

Fig4 <- ggplot(data_line_fin_long, aes(x = median_date, y = support, color = policy_area, group = policy_area)) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = policy_area), alpha = 0.2, color = NA) +
  geom_line(size = 1, na.rm = TRUE) +
  geom_point(size = 2, na.rm = TRUE) +
  labs(x = "", y = "Support (in %)", color = "Policy Area", fill  = "Policy Area") +
  theme_minimal(17) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, size = 16),
    legend.title = element_blank(),
    axis.text.x = element_text(hjust = 1)
  ) +
  scale_x_date(date_labels = "%b %Y", date_breaks = "6 months") +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette  = "Set1") +
  geom_vline(xintercept = as.numeric(as.Date("2022-02-24")), linetype = "dashed", color = "red", size = 1)

png(file = here("Figure4.png"), width = 12, height = 8, units = "in", res = 800)
Fig4
dev.off()

# 8. Figure 6 ############################################################
data_event <- data %>%
  mutate(date = ymd(date, quiet = TRUE))

median_dates <- data_event %>%
  group_by(studyno1) %>%
  summarize(
    median_date = median(date, na.rm = TRUE),
    start_date  = min(date, na.rm = TRUE),
    .groups = "drop"
  )

data_event_fin <- left_join(data_event, median_dates, by = "studyno1") 

data_event_fin <- data_event_fin %>%
  filter(!median_date %in% as.Date(c("2020-11-07", "2021-09-28"))) 
names(data_event_fin)

mod1 <- feols(Proposals_CommonDefence ~  
                i(median_date, closeness_to_russia, ref = as.Date("2022-01-29")) +
                i(median_date, exposure, ref = as.Date("2022-01-29")) +
                gender + age + ageEduc + locality + profession + leftRight + I(leftRight^2)+ 
                polar_eu_position:factor(median_date) +
                weighted_avg_eu_position:factor(median_date) |
                NUTScode + factor(median_date),
              cluster = ~NUTScode , data = data_event_fin)

mod2 <- feols(Proposals_CommonForeign ~  
                i(median_date, closeness_to_russia, ref = as.Date("2022-01-29")) +
                i(median_date, exposure, ref = as.Date("2022-01-29")) +
                gender + age + ageEduc + locality + profession + leftRight + I(leftRight^2)+ 
                polar_eu_position:factor(median_date) +
                weighted_avg_eu_position:factor(median_date) |
                NUTScode + factor(median_date),
              cluster = ~NUTScode, data = data_event_fin)

mod3 <- feols(Proposals_EnlargeEU ~  
                i(median_date, closeness_to_russia, ref = as.Date("2022-01-29")) +
                i(median_date, exposure, ref = as.Date("2022-01-29")) +
                gender + age + ageEduc + locality + profession + leftRight + I(leftRight^2)+
                polar_eu_position:factor(median_date) +
                weighted_avg_eu_position:factor(median_date) |
                NUTScode + factor(median_date),
              cluster = ~NUTScode, data = data_event_fin)

mod4 <- feols(Proposals_CommonEnergy ~  
                i(median_date, closeness_to_russia, ref = as.Date("2022-01-29")) +
                i(median_date, exposure, ref = as.Date("2022-01-29")) +
                gender + age + ageEduc + locality + profession + leftRight + I(leftRight^2)+ 
                polar_eu_position:factor(median_date) +
                weighted_avg_eu_position:factor(median_date) |
                NUTScode + factor(median_date),
              cluster = ~NUTScode, data = data_event_fin)

etable(mod1, mod2, mod3, mod4, 
       tex = TRUE, 
       dict = c(
         "Proposals_CommonDefence" = "Common Defence",
         "Proposals_CommonForeign" = "Common Foreign Policy",
         "Proposals_EnlargeEU" = "EU Enlargement",
         "Proposals_CommonEnergy" = "Common Energy Policy"
       ))

event_plot <- function(model, title_label, y_label) {
  tidy(model, conf.int = TRUE) %>%
    filter(grepl("closeness_to_russia", term)) %>%
    mutate(date = as.Date(gsub(".*::", "", term))) %>%
    bind_rows(
      tibble(
        term = "baseline",
        estimate = 0,
        std.error = NA,
        statistic = NA,
        p.value = NA,
        conf.low = 0,
        conf.high = 0,
        date = as.Date("2022-01-29")
      )
    ) %>%
    ggplot(aes(x = date, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    geom_vline(xintercept = as.Date("2022-01-29"), linetype = "dashed", color = "red", size = 1) +
    geom_pointrange(color = "black") +
    labs(
      x = "",
      y = y_label,
      title = title_label
    ) +
    scale_x_date(
      date_labels = "%b %Y",
      date_breaks = "6 months"
    ) +
    ylim(-0.15, 0.4) + 
    theme_minimal(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_text(angle = 90, hjust = 1, color = "black"),
      axis.text.y = element_text(color = "black"),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
      panel.grid.major = element_line(color = "gray90"),
      panel.grid.minor = element_blank()
    )
}

p1 <- event_plot(mod1, "Defence", "Effect on Support \nfor Common Defence Policy")
p2 <- event_plot(mod2, "Foreign", "Effect on Support \nfor Common Foreign Policy")
p3 <- event_plot(mod3, "Enlargement", "Effect on Support \nfor Enlargement")
p4 <- event_plot(mod4, "Energy", "Effect on Support \nfor Common Energy Policy")

combined_plots <- ggarrange(
  p1, p2, p3, p4,
  ncol = 2, nrow = 2,
  labels = c("A", "B", "C", "D"),
  align = "hv"
)

png(file = here("Figure6.png"), width = 7, height = 7, units = "in", res = 800)
combined_plots
dev.off()

# 9. Figure 7 ############################################################
event_plot <- function(model, title_label, y_label) {
  tidy(model, conf.int = TRUE) %>%
    filter(grepl("exposure", term)) %>%
    mutate(date = as.Date(gsub(".*::", "", term))) %>%
    bind_rows(
      tibble(
        term = "baseline",
        estimate = 0,
        std.error = NA,
        statistic = NA,
        p.value = NA,
        conf.low = 0,
        conf.high = 0,
        date = as.Date("2022-01-29")
      )
    ) %>%
    ggplot(aes(x = date, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    geom_vline(xintercept = as.Date("2022-01-29"), linetype = "dashed", color = "red", size = 1) +
    geom_pointrange(color = "black") +
    labs(
      x = "",
      y = y_label,
      title = title_label
    ) +
    scale_x_date(
      date_labels = "%b %Y",
      date_breaks = "6 months"
    ) +
    ylim(-0.2, 0.25) + 
    theme_minimal(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_text(angle = 90, hjust = 1, color = "black"),
      axis.text.y = element_text(color = "black"),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
      panel.grid.major = element_line(color = "gray90"),
      panel.grid.minor = element_blank()
    )
}

p1 <- event_plot(mod1, "Defence", "Effect on Support \nfor Common Defence Policy")
p2 <- event_plot(mod2, "Foreign", "Effect on Support \nfor Common Foreign Policy")
p3 <- event_plot(mod3, "Enlargement", "Effect on Support \nfor Enlargement")
p4 <- event_plot(mod4, "Energy", "Effect on Support \nfor Common Energy Policy")

combined_plots <- ggarrange(
  p1, p2, p3, p4,
  ncol = 2, nrow = 2,
  labels = c("A", "B", "C", "D"),
  align = "hv"
)

png(file = here("Figure7.png"), width = 7, height = 7, units = "in", res = 800)
combined_plots
dev.off()

